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J> ' We describe the calculation of the two-loop massive operator matrix elements with massive ex- 

ternal fermions in QED. We investigate the factorization of the 0(a 2 ) initial state corrections 
to e + e~ annihilation into a virtual boson for large cms energies s ^> m\ into massive operator 
matrix elements and the massless Wilson coefficients of the Drell-Yan process adapting the color 
coefficients to the case of QED, as proposed by Berends et. al. in Ref. [1]. Our calculations 
show explicitly that the representation proposed in Ref. [ 1 ] works at one-loop order and up to 
terms linear in \n(s/m 2 ) at two-loop order. However, the two-loop constant part contains a few 
structural terms, which have not been obtained in previous direct calculations. 
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Ever since the operator product expansion formalism was applied to the analysis of deep in- 
elastic scattering (DIS), there has been a lot of interest in the calculation of massive operator matrix 
elements at higher order in perturbation theory. Heavy flavor corrections to DIS structure functions 
are very important at small values of the Bjorken variable x (where they contribute on the level 
of 20-40%), and can be calculated in the limit Q 2 > m 2 as a convolution of the corresponding 
massive operator matrix elements and the light flavor Wilson coefficients Here Q 2 denotes the 
virtuality of the gauge boson exchanged in DIS and m is the mass of the heavy quark. 

The scaling violations of the heavy flavor part in the structure functions ^^(x, Q 2 ) are very 
different from those of the light flavor contributions, and their knowledge is very important for 
precision measurements of Aqcd and the extraction of light parton densities. A semi-analytic 
calculation of the heavy flavor contributions was done in Ref. [[|] at next-to-leading order for the 
full kinematic range, and a fast numerical implementation in Mellin space was given in Full 
analytic results in the asymptotic region Q 2 S> m 2 for the structure function F^ix^Q 2 ) at 0(a 2 ) 
were derived in Ref. [^], and recently recalculated in In the same kinematic region, F®®(x, Q 2 ) 
was obtained at 0(q s 3 ) in Ref. [§]. 

The O(a^) massive operator matrix elements are required in both cases and all but 
0((m 2 / Q 2 ) k ), k > 1 contributions are found. In this approximation, the structure function F2(x, Q 2 ) 
turns out to be very well described for Q 2 > 20GeV 2 , while for F L (x, Q 2 ) this approximation only 
holds at large scales Q 2 > lOOOGeV 2 . 

More recently, there has been considerable progress in the calculation of the O(a^) heavy 
flavor contributions to the Wilson coefficients of the structure function F2(x,Q 2 ) and the massive 
gluonic operator matrix elements. First, in Ref. [^] the 0(a 2 e) corrections to these matrix ele- 
ments were given. These corrections are required to perform the corresponding renormalization 
procedure at O(a^). Later, the calculation of these contributions for a number of fixed moments 
N = 2. ..10(12, 14) at O(a^) has also been achieved in [Q], and O(a^rif) contributions were given 
in Ref. [^]. 0(a 2 ) and O(a^) heavy flavor contributions to trans versity have also been obtained in 

On the other hand, massive operator matrix elements, and in particular those with a massive 
external fermion line, can also be applied to a different kind of problem, namely, the calculation of 
initial state QED corrections of scattering processes, such as e + e~ annihilation into a virtual gauge 
boson, using the renormalization group technique. A wealth of information about the Standard 
Model has been obtained in the past from electron-positron colliding beam experiments at different 
facilities around the world. In the future, projects like ILC [11] and CLIC [12] are planned to put 
the Standard Model to even more decisive tests and to reveal new physics [13]. In this context 
high-luminosity machines which operate at a narrow energy regime as DAFNE [ 14 ] and GIGA-Z 
[11, 13] at the Z-peak will offer much higher precision on rare processes. 

The QED initial state radiation causes large corrections for various differential and integral 
scattering cross sections, depending on the sensitivity of the sub-system cross section with respect 
to kinematic rescaling of variables and has to be known at sufficient precision. Both for e + e~ 
annihilation at resonance peaks and the wings of resonances at very high luminosities, the knowl- 
edge of the 0(a 2 ) corrections is mandatory to cope with the experimental precision. While the 
0(a) corrections are known for a large amount of reactions, the corrections beyond the universal 
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contributions 0((aL) k ),\ < k < 5, see [15, 16], to higher orders, were only calculated once at 
two-loop order in Ref. [[]]]. Besides the logarithmic orders 0(a 2 L 2 ,a 2 L) with L = ln(s/m 2 ) and 
m e the electron mass, the constant terms 0(a 2 ) are of interest. 

The renormalization group technique allows to decompose the scattering cross section o(e + e~ 
— > V), V = y*,Z* into massive operator matrix elements and massless Wilson coefficients. The 
fermion mass effects are contained in the former, while the sub-system hard scattering cross sec- 
tions are calculated for massless particles. The corresponding massless Wilson coefficients are 
known from the literature [17, 18] for the Drell-Yan process. In [jl|] this method was used to derive 
all terms up to 0(a 2 L) in addition to the direct calculation. 

The differential scattering cross section can be written in the limit s ^> m 2 as a sum of three 
contributions [j|]: 



dot 



do 



HI 

+ I^JZ + ""e+e- (J) 

ds' ds' ds' ds' ' 

where the labels I, II and III refer to the flavor non-singlet terms with a single fermion line, those 
with an additional closed fermion line, and the pure-singlet terms, respectively. Here, s' denotes 
the invariant mass of the virtual vector boson and s the cms energy of the process. 



da 



■ xs, < x < 1. 



(2) 



It is convenient to write the scattering cross section in Mellin space by applying the integral 
transform 

Using the renormalization group method it can be shown that the three contributions in Eq. ([[]) 
can be expressed as [JTL [l9|] 
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where L = \n(s / m 2 ) +\n{x) , and cr( ) is the Born cross section. The quantities P^ and are the 
LO splitting functions and LO Wilson coefficients, respectively, while at NLO they are denoted by 
PeP and del\ Yj) and F^J are the one-loop and two-loop constant terms of the massive operator 
matrix elements, respectively. The labels I, II and III appear in Pee , deP and r[ l J corresponding 
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Figure 1: The self-energy diagrams. 



to the three possible contributions. The constant j3o is the first coefficient in the expansion of the 
QED j8 -function. 

The splitting functions up to 0{a 2 ) are well known [p0(|, and as we mentioned before, so 
are the massless Wilson coefficients [|l7|, Q. Here we present the missing ingredient needed to 
complete the decomposition of the scattering cross section according to Eqs. (Q-§|), namely, the 
0(a 2 ) massive operator matrix elements. Details on the calculation can be found in Ref. [JlT 

The bare operator matrix elements are given by 



rn- 



-,e,N 



k=i 



,e,N 



(7) 



where a is the unrenormalized coupling constant. The double hat means that the quantity is com- 
pletely unrenormalized. The complete renormalization procedure includes charge renormalization, 
wave function renormalization and the renormalization of the composite operators. We need the 
inclusion of some counterterm diagrams in the case of process I. The electron mass is renormalized 
on-shell, p 2 = m 2 , with p the momentum of the external fermion, which means that no collinear 
singularities will appear. 

The wave function renormalization was performed using the Z-factors coming from the fermion 
self-energy, see Fig. |]. Since we have a massive fermion in the external legs, the coupling constant 
was first obtained in the MOM-scheme, after which we transform to the MS scheme using 
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We perform the calculation in D = 4 + e dimensions. It can be shown that after wave function 
and charge renormalization, keeping the charge in the MOM-scheme, the two-loop OMEs, denoted 
by a single hat, are given by 
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For the final renormalization step, i.e., the renormalization of the composite operators we just 
need to include the corresponding inverse Z-factors, , after which the completely renormalized 
OMEs are given by 
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(f) (g) 00 0) (j) 

Figure 2: Feynman diagrams for the calculation of the massive two-loop operator matrix elements A, 



Figure 3: Counterterm diagrams. The black stars represent the counterterm vertices. 
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The Feynman diagrams required for process I are shown in Fig. ||. As it was mentioned before, 
in order to renormalize the corresponding OMEs for this process, we need to include also the 
counterterm diagrams shown in Fig. |L where the black stars represent the counterterm vertices in 
QED. The diagrams are calculated using the standard Feynman rules for the operator insertions, and 
projecting the resulting numerators in the integrands with the factor (/+m)/4, after which we take 



the trace. The program FORM [21] was used to decompose the diagrams as a linear combination 
of integrals with different powers of propagators, which arise after canceling as many terms as 
possible in the numerator against the propagators. The most complicated integrals appearing in 
the expressions are those where all five propagators are present. These can be computed using 
integration by parts identities to express them in terms of 4-propagator integrals. The resulting 
integrals were checked by several means, including their representation in terms of Mellin-B antes 



integrals, cf. [|19]. 

The constant contributions to process I in Eq. 



is given in x-space by [19] 
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Figure 4: Feynman diagrams for the calculation of the massive two-loop operator matrix elements A 
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where Li„(x) and S„, p (x) are the well known polylogarithm and Nielsen functions, and 
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For process II we need the diagrams in Fig. ||. These are relatively easy to compute, since they 
involve a one-loop insertion. The result for the constant term of process II in Eq. < JTo| ) is given by 
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Figure 5: Feynman diag rams for the calculation of the massive two-loop operator matrix elements A^e ' 



In the case of process III, we need the pure singlet diagrams shown in Fig. |5[ These can be 
calculated using the corresponding one-loop operator insertions, which were given in Ref. [|22|]. 



The result for the constant term appearing in Eq. ( 1 1 ) is [ 19] 
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Our calculations show explicitly that the operator matrix elements satisfy the relations given 
in Eqs. (|9-11) for the pole terms, which automatically guarantees that the decomposition given in 
Eqs. (£§-|6]) holds for the logarithmic terms, as found in Ref. [jl|]. Furthermore, we have verified 
that the first moment of the operator matrix elements vanishes, including the constant terms given 
in Eqs. (13-15|), so they obey fermion number conservation as they should. This provides a highly 
non-trivial check of our results. However, when we assemble the final result for the constant terms 
in Eqs. by including the massless Wilson coefficients, we observe that, although many of 

terms appearing in Ref. ^ appear in our results, a few structural terms, such as terms proportional 
to 1/(1 — x) 3 and ln(x)/(l — x) 4 , also appear, which were not present in the result given in [fi]]. 
This result appears in contrast to the case of massless external fermions and boson lines, where the 
corresponding cross sections have been shown to factorize, including the constant terms, as can 
be seen in Refs. [||, 0, BL 0, 23]. The issue requires further investigation in order to elucidate the 
reasons for this. 
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